r 


AD-A170  051 


SOLUTION  OF  POISSON'S  EOUATION  USING  METHOD  OF  MOMENTS: 
APPLICATION  TO  MOS  DEVICES(U)  RIT  RESEARCH  CORP 
ROCHESTER  NV  E  ARVAS  ET  AL  MAR  86  N00014-84-K-0522 

F/G  20/12 


UNCLASSIFIED 


Ad-  flmo 


!« 

k 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  of  This  PAGE  (When  Data  Entered) 


REPORT  DOCUMENTATION  PAGE 


Wl 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


2.  GOVT  ACCESSION  NO.I  3.  RECIPIENT'S  CATALOG  NUMBER 


.  TITLE  (and  Subtitle) 

SOLUTION  OF  POISSON'S  EQUATION  USING  METHOD  OF 
FOMENTS:  APPLICATION  rC  KOS  DEVICES 


S.  TYPE  OF  REPORT  b  PERIOD  COVERED 

Technical  Resort 
01  Oct.  85  to  31  Mar  86 


6.  PERFORMING  ORG.  REPORT  NUMBER 


AUTHORr*; 


8.  CONTRACT  OR  GRANT  NUMBERS 


N00014-84-K-0532 


.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS  10.  PROGRAM  ELEMENT.  PROJ  ECT,  T  A5K 

AREA  ft  WORK  UNIT  NUMBERS 

RIT  Research  Corporation 

75  Highpower  Road  NR  613-005 

Rnrhpster.  NY  14623-3435 _ 

t.  CONTROLLING  OFFICE  NAME  AND  ADDRESS  12.  REPORT  DATE 

Department  of  Navy  March  1986 

Office  of  Naval  Research  13.  NUMBER  OF  PAGES 

Arlington,  VA  22217 _ _ _ 33 _ 

l«.  MONITORING  AGENCY  NAME  «  ADDRESSf/I  dtllerant  trom  Controlling  Olttce)  IS.  SECURITY  CLASS,  (ol  tbit  report) 

UNCLASSIFIED 


Ercument  Arvas,  Perambur  S.  Neelakantaswamy,  and 
Ibrahim  R.  Turkman 


NR  613-005 

12.  REPORT  DATE 


13.  NUMBER  OF  PAGES 

33 


1S«.  DECL  ASSIFICATION/OOWNCRADING- 
SCHEDULE 


j _ 

16  niCTRiBiiTION  STATEMENT  (ot  thle  Report) 


Thin  cl.  _ua:»v„ 

for  pvbU-t  r*Jiwr-.  cc*i  «ot«  U 
dittiribv’ioa  fa  4 


M7.  DISTRIBUTION  STATEMENT  (of  tn»  aoxrcu  «*«.« 


It  different  trom  Report) 


jUL  0  7  W86 


118-  SUPPLEMENTARY  NOTES 


19.  KEY  WORDS  (Continue  on  ravaree  aide  II  neceeemry  mid  Identify  by  block  number)  | 

M0S  DEVICES,  POISSON'S  EQUATION  (TWO-DIMENSIONAL),  METHOD  OF  MOMENTS, 

THRESHOLD  VOLTAGE,  SURFACE  POTENTIAL 

20.  ABSTRACT  fCondnu.  on  reverme  etoe  It  neceeemry  end  Identllv  by  block  number)  . 

An  algorithm  for  the  computation  of  solution  to  Poisson's  equation  in  a  two- 
dimensional  domain  is  developed  in  terms  of  equivalent  sources  on  the  boundary 
The  region  considered  can  be  of  arbitrary  shape,  and  the  boundary  conditions 
can  be  Dirichlet,  Neumann  or  mixed  type.  The  solution  is  obtained  by  method 
of  moments.  Pulse  expansion  and  point  matching  techniques  are  used.  Computed 
results  closely  agree  with  the  available  data  concerning  M0S  devices. 


DD  ,  :°r7,  1473  EDITION  OF  I  NOV  SS  IS  OBSOLETE 


S/N  0102-LF-014-6601 


_ UNCLASSIFIED _ 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  flTh.n  Data  Entered) 


i  - .  .  - 


SOLUTION  OF  POISSON'S  EQUATION  USING  METHOD 
OF  MOMENTS:  APPLICATION  TO  MOS  DEVICES 


by 


Ercument  Arvas 

Electrical  Engineering  Department 
Rochester  Institute  of  Technology 
One  Lomb  Memorial  Drive 
Rochester,  New  York  14623 


Perambur  S.  Neelakantaswamy 
RIT  Research  Corporation 
75  Highpower  Road 
Rochester,  New  York  14623 


and 


Ibrahim  R.  Turkman 
Electrical  Engineering  Department 
Rochester  Institute  of  Technology 
Rochester,  New  York  14623 


RIT  Research  Corporation 
75  Highpower  Road 
Rochester,  New  York  14623 


March,  1986 


l  i:.:»  . 

4  tor 

I  viUtrlbv’UH. 


1 


CONTENTS 

Page 

ABSTRACT  .  1 

I.  INTRODUCTION .  2 

II.  STATEMENT .  3 

III.  METHOD  OF  SOLUTION .  3 

IV.  SAMPLE  RESULTS .  4 

V.  APPLICATION  TO  MOS  STRUCTURES  .  5 

VI .  REFERENCES .  6 

APPENDIX  A:  COMPUTER  PROGRAM .  8 

PROGRAM  LISTING .  11 

SAMPLE  INPUT/OUTPUT  FILE  .  23 

FIGURES  1-8 .  26 


ABSTRACT 


An  algorithm  for  the  computation  of  solution  to  Poisson's 
equation  in  a  two-dimensional  domain  is  developed  in  terms  of 
equivalent  sources  on  the  boundary.  The  region  considered 
can  be  of  arbitrary  shape,  and  the  boundary  conditions  can  be 
Dirichlet,  Neumann  or  mixed  type.  The  solution  is  obtained 
by  method  of  moments.  Pulse  expansion  and  point  matching 
techniques  are  used.  Computed  results  closely  agree  with  the 
available  data  concerning  MOS  devices. 


I. 


INTRODUCTION 


Poisson's  equation  is  one  of  the  most  important 
differential  equations  of  physics.  For  example,  it  can  be 
used  to  find  the  threshold  voltages  of  MOSFET's.  When  the 
channel  length  is  small,  the  depletion-layer  widths  of  the 
source  and  drain  junctions  are  comparable  to  the  channel 
length,  and  the  potential  distribution  is  two  dimensional 
amenable  for  solution  via  Poisson's  equation. 

In  this  work  we  give  a  simple  method  for  solving  two- 
dimensional  Poisson's  equation  in  a  region  subject  to  general 
boundary  conditions  on  the  bounding  surface.  Equivalent 
surface  charges  are  placed  just  outside  the  boundary  and  the 
total  potential  (produced  by  the  impressed  volume  charges  and 
the  equivalent  surface  charges)  is  enforced  to  satisfy  the 
boundary  conditions.  This  transforms  the  boundary  value 
problem  into  an  integral  equation  for  the  equivalent  surface 
charges.  Then  the  method  of  moments  [1]  is  used  to 
solve  the  integral  equation  numerically. 


II. 


STATEMENT  OF  PROBLEM 


Consider  a  2-dimensional  region  R  bounded  by  the  contour 
C  as  shown  in  Figure  1.  The  problem  is  to  find  the  total 
potential  ip  (x,y)  in  R  which  satisfies  the  Poisson's  equation 


9  2 1 p  +  3  2 

3x  2  3y  2 


=  -pv/e 


(1) 


in  R,  with  the  boundary  condition(s) 

+  8  H  =  y 


(2) 


on  C. 

In  eqn.  (1),  Pv  denote  the  volume  charge  density,  and  e 
is  the  permittivity  of  the  material  in  R. 

I 

In  eqn.  (2),  a  ,  B  and  y  denotes  known  functions  defined 

on  C. 


Note  that  the  general  condition  of  eqn.  (2)  includes,  as 
special  cases,  of  Dirichlet  (a  -1,  B  =o)  and  Neumann  (  a-o,  B 
-1)  conditions. 

The  Laplace's  equation  is  the  special  case  of  eqn.  (1) 
with  p  -  o  The  solution  of  Laplace's  equation  is  given  in 
detail  ii7  f 2 ] .  The  present  work  is  an  extension  of  the  work 
in  [2],  modified  for  Poisson's  equation  applied  to  MOS 
structures . 


III.  METHOD  OF  SOLUTION 


In  solving  eqn.  (1)  subject  to  boundary  condition  of 
eqn.  (2)  we  let 


+  *h 


(3) 


where 


(in  R) 


(4) 


V  ‘  \  =  0  (in  R) 


and 


(5) 


k 

r  -  r' 


dx'dy’ 


(6) 


The  solution  to  eqn.  (4)  is 


<P 


P 


_ n 

2t:£  r 


p(r')  In  — 


where  £  and  £'  denote  the  radius  vector  to  the  field  and  a 
source  point  respectively,  p  (£' )  is  the  value  of  the 
impressed  charged  density  at  r's  and  k  is  an  arbitrary 
constant  (taken  as  100.0  in  tKis  work). 

The  Laplacian  potential  can  be  assumed  to  be  produced 
by  some  equivalent  surface  charges,  a  ,  outside  R  (Fig.  2).. 
Hence  j 0^  is  the  solution  of  eqn.  (5)  subject  to  boundary 
condition 

34>h  (7) 

+  e  =  Y'a^p  "  e  a/  on  c 

Since  £Th  has  the  form 

'  1  (8) 

+h  -  1¥T  c  °  ln  ~r  -  r'  ix'dy' 


we  see  that  eqn.  (7)  is  an  integral  equation  for  a  . 

Note  that  (5)  subject  to  the  boundary  condition  of  eqn. 

(  7)  is  the  same  boundary  value  problem  as  the  one  considered 
in  [2].  We  use  pulse  expansion  and  point-matching  techniques 
to  solve  this  problem. 

+ 

The  approach  involved  is  to  first  model  the  surface  C 
by  N  planar  strips  and  the  assume  a  constant  charge  density 
on  each  segment.  Satisfying  the  boundary  condition  of  eqn. 

(7  )  at  the  center  of  each  of  N  strips,  gives  N  algebraic 
equations.  The  solution  of  these  equations  gives  the  value 
of  the  constant  charge  density  on  each  strip.  The  details 
are  elaborated  in  (2).  Once  eqn.  (7)  is  solved  for  o  ,  we 
obtain  the  total  potential  \p  using  eqns.  (8),  (6)  and  (3). 


IV.  SAMPLE  RESULTS 


A  FORTRAN  program  is  written  to  implement  the  theory 
developed  above. 


The  first  test  problem  that  we  tried  is  shown  in  Fig. 

-12 

3,  where  a  line  charge  of  p  -  8.854x10  C/m  is  placed  at 

1 

the  center  of  a  grounded  rectangular  boundary.  The  total 
potential  was  evaluated  at  the  points  A,  B,  C  and  D  as  shown. 

Table  1  illustrates  the  convergence  of  the  computed 
results  as  the  number  N  of  segments  is  increased.  The  exact 
result  [3,  eqn.  4-7.23]  is  also  shown  for  comparison.  The 
last  column  of  the  table  shows  the  CPU  time  on  a  VAX  11/782. 

A  second  test  problem  formulated  to  study  a  short 
channel  MOSFET  is  described  hereunder. 


VI.  APPLICATION  TO  MOS  STRUCTURES 

To  demonstrate  the  applicability  of  the  proposed 
numerical  method  to  MOS  structures,  a  test  N-channel  MOSFET 
illustrated  in  Fig.  4  is  considered.  The  rectangular 
depletion  region  under  the  gate  and  its  expanded  view  with 
the  relevant  boundary  conditions  are  depicted  in  Fig.  5. 

The  notations  followed  are  those  detailed  in  [4],  Figs. 
6  and  7  illustrate  the  surface  potential'!'  (x,  y  -  d) 
variation  along  the  channel  for  2  typical  devices  with 
channel  lengths  L  ■  lu  m  and  5  y  m  respectively. 

The  corresponding  threshold  voltage  (V  -v  )  versus 

T  FB 

channel  length  for  drain  voltage  (V  )  of  0  and  5V  is 
presented  in  Fig.  8.  D 

For  comparison,  along  with  the  computed  data,  the 
results  obtained  by  (approximate)  closed-form  solution  due  to 
Poole  and  Kwong  [4]  are  also  shown  in  Figs.  6,  7  and  8. 

Referring  to  these  figures  (Figs.  6,  7  and  8),  close 
agreement  between  the  results  may  be  observed.  Any  deviation 
can  be  attributed  to  the  approximations  involved  in  the 
truncation  of  the  series  solution  given  in  (4)  and  due  to  the 
variations  in  the  values  of  d  and  V  considered  in  the 
analysis.  gm 

However,  the  present  work  indicates  the  applicability  of 
the  method  of  solution  envisaged  to  the  MOS  structures.  This 
method  can  be  extended  to  a  more  realistic  model  of  the  MOS 
structure  involving  curved  depletion  boundaries  and  the 
depletion  width  (d)  varying  along  the  channel  length. 


Further,  this  steady  state  solution  can  be  extended  to 
study  transient  causes  pertaining  to  ESD/EOS  induced  effects. 
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Table  I 


Potential  computed  for  the  problem  of  Fig.  3 


N 

A 

(0.25,0.75) 

B 

(0.25,0.50) 

C 

(0.1,0.25) 

D 

(0.45,0.95) 

CPU  Time 
( Sec . ) 

4 

0.05924 

0.1133 

0.01411 

0.3145 

3.28 

8 

0.07271 

0.1268 

0.02757 

0.3280 

3.55 

12 

0.07113 

0.1226 

0.02875 

0.3244 

4.03 

16 

0.07050 

0.1221 

0.02747 

0.3238 

4.54 

20 

0.07032 

0.1219 

0.02721 

0.0236 

5.15 

24 

0.07024 

0.1218 

0.02719 

0.3235 

5.92 

32 

0.07018 

0.1217 

0.02714 

0.3234 

7.86 

40 

0.07016 

0.1217 

0.02712 

0.3234 

10.69 

60 

0.07014 

0.1216 

0.02711 

0.3234 

21.43 

80 

Exact 

0.07014 

0.1216 

0.02711 

0.3234 

37.59 

APPENDIX  A:  COMPUTER  PROGRAM 


The  FORTRAN  computer  program  is  composed  of  a  main  and  9 
subprograms.  The  subprograms  are: 

INFOR 

SOLTN 

VMATRX 

ZMATRX 

FIELD 

ELSV 

POTEN 

INTG 

GRAD 

The  last  three  programs  compute  the  potential  and  its 
gradient  at  a  point  (x,y)  due  to  the  impressed  charge 
distribution.  Hence,  as  the  source  is  charged,  these 
programs  must  be  changed  accordingly. 

The  Main  Program: 

The  main  program  reads  in: 

a)  the  number  ( NTOTAL )  of  the  straight  line  segments 
approximately  the  boundary  C  of  the  region  R. 

b)  the  dielectric  constant  (SPSR)  of  the  medium  R. 

c)  the  parameter  LAPOIS.  If  LAPOIS  is  equal  to  zero,  the 
problem  is  to  solve  the  Laplace's  equation.  (In  this 
case  the  last  three  subroutines  are  not  needed).  If 
LAPOIS  is  equal  to  1,  we  are  solving  Poisson's  equation 
and  hence  the  potential  and  its  gradient  produced  by  the 
impressed  sources  must  be  provided  by  the  last  three 
subroutines. 

For  each  of  NTOTAL  linear  segments,  the  main  program 
calls  the  subprogram  INFOR.  Then  it  calls  the  subroutine 
SOLTN. 

The  INFOR  subprogram: 

The  subroutine  INFOR  reads  in 

a)  The  coordinates  (Xl,  Yl)  and  X2,  Y2)  of  the  starting  and 
ending  points  of  each  linear  segment  approximating  the 
boundary  C.  (Xl,  Yl,  X2,  Y2  are  in  micrometers). 

b)  The  number  NSEC,  of  subsections  that  each  linear  segment 
is  to  be  divided  into. 


c)  For  each  linear  segment  a  ,  3  and  a  are  read  in.  These 
are  sent  back  to  the  main  program,  where  they  are  stored 
in  the  matrix  BCOND. 

In  the  subroutine  INFOR,  the  coordinates  of  the  starting 
and  ending  points  of  each  subsection  is  computed.  This 
information  is  stored  in  arrays  XVl ,  YVl,  XV2  and 

The  subroutine  SOLTN : 


In  this  subroutine  the  moment  matrix  equation  is  formed 
and  solved.  This  subroutine  calls  various  subroutines. 

i )  The  subroutine  VMATRX : 

In  this  subroutine  the  right  hand  side  of  eqn.  (7)  is 
computed  at  the  center  of  each  subsection.  The  result  is 
stored  in  the  array  V. 

ii )  The  subroutine  ZMATRX : 

In  this  subroutine  the  moment  matrix  Z  is  computed.  The 
(i,j  )  th  element  of  this  matrix  is  the  right  hand  side  of 
eqn.  (7),  computed  at  the  center  of  jth  subsection.  (0  here 
is  the  potential  produced  by  a  constant  charge  density  of 
2 

l(C/m  )  on  the  jth  subsection). 

i i i )  The  subroutine  ELSV: 

This  subroutine  takes  the  inverse  of  the  moment  matrix  Z 
and  stores  the  inverse  matrix  into  the  Z  matrix. 

Once  Z  matrix  is  inverted,  the  surface  density  is 
computed  in  SOLTN  subroutine  by  multiplying  the  inverse  of 
the  Z  matrix  by  tne  column  vector  V.  The  charge  density  is 
stored  in  the  array  I. 

iv)  The  subroutine  FIELD : 


This  subroutine  computes  the  total  potential  at  K  points 
equally  spaced  between  the  points.  (XIN,  YIN)  and  XFIN,  YFIN ) 

The  last  three  subroutines  compute  the  potential  and  its 
gradient  at  a  given  point  due  to  a  constant  volume 

3 

charge  density  RHO  (C/m  )  in  a  rectangular  cylinder  of 
infinite  length. 
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INPUT/OUTPUT  OF  THE  PROGRAM 


The  input  to  the  program  is  through  the  data  file  92. 

The  first  line  of  the  input  file  is 

NTOTAL ,  EPSE,  LAPOIS 

Then  we  have  NTOTAL  pairs  of  lines  which  have  the  form 
XI,  Yl,  X2,  Y2,  NSEC 

ALPHA,  BETA,  GAMA. 

Where  (XI,  Yl)  and  (X2,  Y2 )  denote  the  coordinates  of  the 
starting  and  ending  point  of  a  linear  segment,  and  NSEC  is 
the  number  of  subsection,  that  the  segment  will  be  subdivided 
into  ALPHA,  BETA  AND  GAMA  show  the  values  of  a  ,  B 
and  o  on  the  segment. 

The  last  line  in  the  input  file  has  the  form 
XIN,  YIN,  XFIN ,  YFIN ,  K  where  (X  IN,  YIN)  and  (XFIN,  YFIN ) 
are  the  coordinates  of  two  points,  and  K  is  an  integer.  The 
program  will  compute  the  total  potential  at  K  equidistant 
points  lying  between  the  points  (XIN,  YIN)  and  (XFIN,  YFIN). 

The  output  of  the  program  is  printed  in  data  file  18. 
Here  the  potential  at  K  points  is  printed.  X  and  Y  are  the 
coordinates  of  the  point  at  which  the  potential  is  computed. 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


C  THIS  PROGRAM  COMPUTES  THE  EQUIVALENT  ELECTRIC  C 
C  CHARGE  DENSITY  ON  THE  SURFACE  OF  A  LOSSLESS  DIELECTRIC  C 
C  CYLINDER.  THIS  IS  A  TWO-DIMENSIONAL  PROBLEM.  THE  TOTAL  C 
C  POTENTIAL  INSIDE  THE  DIELECTRIC  IS  DUE  TO  SOME  SPECIFIED  C 
C  CHARGES  INSIDE  AND  DUE  TO  SOME  IMPRESSED  POTENTIALS  C 
C  ALONG  THE  SURFACE.  C 
C  C 
C  AT  ANY  POINT  ON  THE  BOUNDARY  OF  THE  CYLINDER  WE  HAVE  C 
C  C 
C  ALPA( C ) *POT( CHARGE ) +BETA{ C ) * { D/DN ) ( POT ( CHARGE ) )  =  -ALPHA(C)  C 
C  *POT ( SOURCE ) -BETA( C )  *  ( D/DN ) POT ( SOURCE ) +GAMA( C )  C 

c  c 

C  WHERE;  C 
C  C  SHOWS  THE  VARIABLE  ALONG  THE  BOUNDARY  OF  THE  CYLINDER,  C 
C  ALPHA ( C ) ,  BETA(C)  AND  GAMA(C)  ARE  THREE  FUNCTIONS  THAT  ARE  C 
C  SPECIFIED  AT  ANY  POINT  C,  C 
C  POT ( CHARGE )=POTENTIAL  PRODUCED  AT  THE  POINT  C,  BY  THE  C 
C  UNKOWN  EQUIVALENT  SURFACE  CHARGE  RESIDING  ON  THE  BOUNDARY  C 
C  OF  THE  CYLINDER,  C 
C  (D/DN)  IS  AN  OPERATOR  WHICH  GIVES  THE  NORMAL  DERIVATIVE  C 
C  OF  THE  FUNCTION  THAT  IT  OPERATES  ON  ,  AND  C 
C  POT( SOURCE)  IS  THE  POTENTIAL  PRODUCED  BY  THE  IMPRESSED  C 
C  SOURCES  AT  THE  POINT  ON  THE  BOUNDARY.  THE  IMPRESSED  C 
C  SOURCES  ARE  THE  VOLUME  CHARGE  DENSITY  INSIDE  THE  CYLINDER.  C 
C  THESE  ARE  THE  SOURCES  THAT  APPEAR  ON  THE  RIGHT  C 
C  HAND  SIDE  OF  THE  POISSON'S  EQUATION.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

C 

IMPLICIT  COMPLEX*16  (C) 

IMPLICIT  REAL*8 ( A-B , E-H , P-Z ) 

DIMENSION  V(200) ,RI(200) ,Z( 20 0,200) 

DIMENSION  XVI (200) ,YV1(200) ,XV2(200) ,YV2(200) 

COMMON/NNN/NDP ( 10 ) ,BCOND(10,3) , NTOTAL 
COMMON/TYPE/LAPOIS , EPSR 
NMAX=200 
C 

C  READ  TOTAL  NUMBER  , NTOTAL, OF  LINEAR  SEGMENTS  WHICH 
C  CONSTITUTE  THE  BOUNDARY  FOR  THE  PROBLEM.  ALSO 
C  READ  THE  DIELECTRIC  CONSTANT  ,EPSR,  OF 
C  THE  MEDIUM. 

C 

C  IF  LAPOIS  IS  ZERO  THEN  WE  ARE  SOLVINC  LAPLACE  EQUATION 
C  IF  LAPOIS  IS  ONE  THEN  THE  PROBLEM  IS  POISSON  TYPE. 

C 

READ (92,*)  NTOTAL, EPSR, LAPOIS 
I F ( LAPOIS . EQ . I ) WR ITE (  9  3 , 1 2  3  2  ) 

IF ( LAPOIS. EO. 0) WRITE ( 93, 1233) 

1232  FORMAT(/,25X, 'THIS  IS  POISSON  S  EQUATION:',/) 

1233  FORMAT(/,25X, 'THIS  IS  LAPLACE  S  EQUATION:’,/) 

WRITE( 93 , 1234 ) NTOTAL, EPSR 

1234  FORMAT(  2  5X,  ’ - ' 

4  ,// , 10X , ' NO .  OF  TOTAL  LINEAR  SEGMENT  BOUNDARI ES= '  ,  1 2  ,/ , 1 OX , 

4  'THE  DIELECTRIC  CONSTANT  OF  THE  CYLINDER  IS= ' , FB . 4 ,//) 

N  A I  =  0 
NB  1  =  0 
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c 

DO  199  1  =  1,  NTOTAL 
C 

C  FOR  EACH  OF  NTOTAL  LINEAR  SEGMENTS  FORMING  THE  BOUNDARY 
C  CALL  THE  INFORMATION  (INFOR)  SUBROUTINE  TO  ; 

C  A)  READ  IN  THE  COORDINATES  ( X 1 , Y 1 )  OF  THE  INITIAL  POINT  AND  (X2,Y2) 

C  OF  THE  FINAL  POINT  OF  I  ’  TH  LINEAR  SEGMENT. 

C  B)  READ  IN  THE  NUMBER  ,NSEC,  OF  SMALLER  SUBSECTIONS  THAT  THIS  PARTICULAR 
C  LINEAR  SEGMENT  IS  TO  BE  DIVIDED.  THE  CENTER  OF  EACH  OF  THESE 
C  SUBSECTIONS  IS  A  MATCHING  POINT. 

C  C)  READ  THE  VALUES  ALPHA,  BETA,  AND  GAMA  FOR  THIS  PARTICULAR 
C  LINEAR  SEGMFNT . 

C  D)  FIND  THE  COORDINATES  OF  THE  STARTING  AND  ENDING  POINTS  OF  THESE 
C  SUBSECTIONS  AND  STORE  THEM  IN  THE  ARRAYS  XV 1 , YVl ,XV2 , YV2 . 

C 

c 

c 

CALL  INFOR (XVI , YVl , XV2 , YV2 , MAI , UMAX , A , B ,C ) 

NDP ( I ) =NAI-NB I 
BCCND ( I , 1 )=A 
BCOND ( I , 2 ) =B 
BCOND ( I , 3 ) =C 
C 

C  WRITE  THE  BOUNDARY  CONDITIONS  DATA  FOR  THIS  LINEAR  SECMENT; 

C 

WRITE (93,111) 

111  FORMAT ( ' 1 1 ) 

WRITE (93,112)  I , A , B , G 

112  FORMAT (////5X, 'THIS  IS  THE  INFORMATION  OF  THE 
$  BOUNDARY  = ' , IX, 13 ,// , 5X, ' HERE  ALPHA3 ' , F9 . 5 , 3X , 

$  'BETA3 ' ,E11.4,3X, 'GAMA= ' ,F9.5,/) 

C 

C  write  GEOMETRICAL  DATA  FOR  THIS  LINEAR  SEGMENT; 

C 

WRITE (93,114) 

114  FORMAT ( / ///l 2X , 2 ( ' X-COORD INATE ’ , 5X , ’ Y-COORDINATE ’ , 5X ) ) 

WRITE ( 93 , 115) (J,XV1(J) , YV 1 ( J ) , XV2 ( J ) ,YV2(J) ,J=NBI+1 ,NAI) 

115  FORMAT ( //( 5X , 1 3 , 4 ( 2X , IE ) ) / ) 

NBI=NAI 

199  CONTINUE 

C 

C  OBTAIN  THE  TOTAL  NUMBER  OF  UNKNOWNS  IN  THE  MATRIX  EOUATION. 

C 

NUNKNS=NAI 

WRITE  (93,993  )  NUNKNS 

993  FORMAT( 5X, 'TOTAL  NO.  OF  UNKNOWNS3 ',13 ,// ) 

C 

C  CALL  THE  SOLUTION  SUBROUTINE  TO  SOLVE  THE  PROBLEM. 

C 

CALL  SOLTN (  Z  ,  V  ,  P.I  ,  XVI  ,  YVl , XV2 , YV2 , NUNKNS , NMAX ) 

998  CONTINUE 

STOP 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  INFOR ( XVI , YVl ,XV2 , YV2 , NAI ,NM,A,P ,0) 

IMPLICIT  COMPLEX* 16  (C) 

IMPLICIT  PEAL*8 ( A-B , E-H , P-Z ) 

C 

C  IN  THIS  SUBROUTINE  THE  DATA  IS  ARPANCED  IN  THE  PROPER  FCPM 
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C  FOP  FURTHFR  COMPUTATIONS. 

C 

DIMENSION  XVI (NM) ,  YVl(NM) ,XV2(NM) ,YV2(MM) 

NNODES=HAI 

C  READ  THE  COORDINATES  (Xl,Yl)  AND  (X2,Y2), 

C  P.EAD  THE  NUMRER  OF  SECTIONS  OF  THE  BOUNDARY  (NSEC), 

C  ALSO  READ  THE  BOUNDARY  CONDITIONS  INFORMATION;  ALPHA(A) , 

C  BETA ( B )  AND  GAMA(C) 

C 

C 

R£AD( 9  2 , * ) X 1 , Yl ,X2,Y2,NSEC 

XI =X1 *  1 . D-0  6 

Y1=Y1*1 .D-06 

X2=X2*1 . D-06 

Y2= Y2* 1 . D-0  6 

READ( 92 ,* )A,B,G 

EDELX=( X2-X1 ) /FLOAT ( NSEC ) 

EDELY=(Y2-Y1 ) /FLOAT ( NSEC ) 

DO  20  J= I , NSEC 
NNODES=NNODES+l 

XVI ( NNODES ) =X 1+FLOAT (J-l ) *  EDELX 
YV 1 ( NNODES ) =Y 1+FLOAT (J-l) *  EDELY 
20  CONTINUE 

DO  70  I=NAI+1 ,NNODES-l 
XV2 ( I) =XV1 ( 1  +  1 ) 

YV2 ( I ) =  YVl (1+1) 

70  CONTINUE 

75  XV2( NNODES )=X2 
YV2( NNODES )=Y2 

76  NAI=NNODES 
RETURN 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  SOLTN ( 2 , V , RI , XV 1 , YV1 , XV2 , YV2 , N , NM ) 

C 

C  IN  THIS  SUBROUTINE  THE  MATRIX  EQUATION  AX=Y  IS  SOLVED  USING  THE 
C  METHOD  OF  MOMENTS. 

C 

IMPLICIT  COMPLEX* 1 6  (C) 

IMPLICIT  REAL*8(A-B,E-H,P-Z) 

DIMENSION  V(N)  ,RI(N)  ,Z(N,N)  ,AUXl(600)  ,AUX2(600) 

DIMENSION  XVI (NM) ,YV1(NM) ,XV2(NM) ,YV2(NM) 

COMMON/NNN/NDP (10) ,BCOND(10,3) , NTOTAL 
C 

C  INTIALIZE  THE  VECTORS  Z,V,AND  RI . 

C 

DO  5  I  =  1 ,  N 
V( I ) =0 . DO 
R I  (  I ) =0 . DO 
5  CONTINUE 

DO  10  1=1, N 
DO  10  J  =  1  ,  N 
Z(  I  ,  J  )  =0  .  DO 
10  CONTINUE 
C 

C  CALL  THE  SUBROUTINE  VMATRX  TO  COMPUTE  THE  EXCITATION  VECTOR. 

C 

CALL  VMATPX ( V , N , XV 1 , YV1 ,XV2 , YV2 ,NM) 

C 
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C  CALL  THE  7.MATRX  SUBROUTINE  TO  OETAIN  THE  IMPEDANCE  MATRIX. 

C 

CALL  ZMATRX ( Z,XV1 , YV 1 , XV2 , YV2 , M , MM ) 

C 

C  CALL  THE  ELSV  SUBROUTINE  TO  INVERT  THE  MATRIX. 

C 

EP=0 . ID-09 

CALL  ELSV{ Z , AUXl , AUX2 ,M  , DE , EP ) 

WRITE ( 93 , 1 18 ) PE 
118  FORMAT ( 5X , ' DE  =’,1E) 

C 

C  MULTIPLY  THE  INVERSE  OF  Z-MATRIX  WITH  THE  EXCITATION  VECTOR 
C  TO  OBTAIN  THE  CHARGES. 

C 

DO  25  1=1, N 
SUM=0 . DO 
DO  24  J  = 1 , N 
SUM=SUM+Z< I,J)*V(J) 

24  CONTINUE 
RI( I)=SUM 

25  CONTINUE 
C 

C  WRITE  THE  CHARGES  ON  THE  OUTPUT  FILE. 

C 

NF=0 

DO  135  JKLM=1 , NTOTAL 
NI=NF+1 

MF=MF+NDP ( JKLM ) 

WRITE( 93,101) 

101  FORMAT ( ' 1 ' ) 

WRITE (93 ,102)  JKLM 

102  FORMAT (//5X, 'CHARGES  ON  THE  BOUND.  =’,1X,I2,//) 

ALPHA=BCOND ( JKLM , 1 ) 

BETA=BCOND( JKLM, 2) 

GAMA=BCOND( JKLM , 3 ) 

WRITE (9 3, 198 ) ALPHA , BETA , GAMA 
198  FORMAT (5X, 'FOR  THIS  BOUNDARY  ALPHA= ' , El 1 . 4 , 

&  3X , ' BETA= '  , El 1 . 4 , 3X  ,  ’ GAMA= ’  , El 1.4,//) 

DO  30  I=NI,NF 
WRITE (93,105)  I , RI ( I ) 

105  FORMAT ( /IX, 13, 5X, Ell. 4) 

30  CONTINUE 
135  CONTINUE 

CALL  FIELD{ RI ,N,XV1 ,XV2 , YV1 , YV2 ,NM) 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  VMATRX ( V , N ,XV1 , YV1 ,XV2 ,YV2,NM) 

C 

C  IN  THIS  SUBROUTINE  THE  EXCITATION  VECTOR  V  IS 
C  COMPUTED  . 

C 

IMPLICIT  COMPLEX* 1 6  (C) 

IMPLICIT  REAL* 8 ( A-B , E-H , P-Z ) 

DIMENSION  V ( N )  ,XV1(NM)  ,YVl(NM)  ,XV2(NM)  ,YV2(NM) 
COMMON/NNN/NDP ( 10) ,BCOND( 10,3) , NTOTAL 
COMMON/TYPE/LAPOIS , E PSP 
PI=4 . DO* DAT AN ( 1 .DO ) 

TP=2 . D0*PI 
EPS=EPSP*8 .854D-12 
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non  .  nnn 


if(lapois.eo.1)tp=tp*fps 

NF  =  0 

DO  lin  JKLMN=1,NT0TAL 
NI=NF+1 

NF“NF+NDP (JKLMN ) 

ALPHA*BCOND  ( JKLMN  ,  1 ) 

8ETA=*BC0ND(  JKLMN,  2) 

GAMA  =  BCOND( JKLMN  ,3 ) 

DO  100  I=N  I , NF 
V(  I ) =GAMA*TP 
IF ( LAPOIS . EQ . 0 ) CO  TO  100 
X1=XV1( I) 

Yl* YV1 ( I ) 

X2=XV2( I) 

Y2=YV2( I) 

XF= ( X1+X2 ) /2 . 0 
YF={ Y1+Y2J/2.0 
IF(ALPHA.EQ.O.O)GO  TO  90 
CALL  POTEN ( XF , YF , POT ) 

V( I ) =V ( I ) -ALPHA* POT 

C  IF(  LAPOIS  .  EQ  .  1  .AND  .BETA .  EO  .0.0)V(I)=*10. 0*V  (  I ) 

IF(BETA.£Q. 0.0)CO  TO  100 
90  CALL  GRAD( XF , YF , POTX , POTY ) 

FNRMD=- ( X2-X1 ) *  POTY 

FNRMD*FNRMD+( Y2-Y1 )*POTX 

FL=SQRT( (X2-X1)*(X2-X1)+(Y2-Y1)*(Y2-Y1) ) 

V( I ) =V ( I ) -BETA*FMRMD/FL 
100  CONTINUE 
110  CONTINUE 

RETURN 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  2MATRX( Z ,XV1 , YV1 ,XV2 , YV2 , N , NM ) 

IN  THIS  SUBROUTINE  THE  Z-MATRIX  IS  FORMED. 

IMPLICIT  COMPLEX* 16  (C) 

IMPLICIT  REAL*8(A-B,E-H,P-Z) 

DIMENSION  XVl(NM) ,YV1(NM) ,XV2(NM) ,YV2(NM) ,Z(N,N) 
COMMON/NNN/MDP ( 10) ,BCOND(10,3) , NTOTAL 
COMMON /TYPE/ LAPOIS , FPSR 
Cl-DCMPLX ( 1 . DO  f  0 . DO ) 

CK=DCMPLX( 100. DO, 0. DO) 

PI=4 . D0*DATAN ( 1 . DO ) 

NF=0 

DO  1000  JKLM=1, NTOTAL 
ALPHA=BCOND ( JKLM , 1 ) 

BETA=BCOND( JKLM , 2 ) 

NI=NF+1 

NF=NF+NDP( JKLM) 

DO  999  I»NI,NF 

COMPUTE  THE  PARAMETERS  OF  THE  FIELD  SUBSECTION 

XI=XV1( I) 

XIP1=XV2( I) 

YI»YV1( I) 

YIP1=YV2( I ) 

CZ I=CMPLX ( X I , YI ) 

CZ IP  1  =  DCMPLX ( X I P 1 ,YIP1 ) 


EDFLI=CDA3S(CZIP1-CZI) 

CUI=(CZIP1-CZI)/EDELI 
CZICAR=( CZI+CZIP1 )/2 -  DO 
DO  888  J= 1 , N 
C 

C  COMPUTE  THE  PARAMETERS  OF  THE  SOURCE  SUBSECTION 
C 

xj»xvi( J) 

XJP 1*XV2 ( J ) 

Y  J  =  YV 1  ( J  ) 

YJP1=YV2 ( J ) 

CZJ*DCMPLX(XJ , YJ ) 

CZJP1 =DCMPLX ( XJPI , YJP1 ) 

EDELJ =CDABS (CZJPL-CZJ) 

CUJ  ={C2JP1-CZJ ) /EDELJ 

CARC= ( CZICAR-CZ JPI ) / ( CZ ICAR-C2 J ) 

C 

CTERM=CDLOG ( CARG ) /CUJ 

IF ( BETA .EO.O.DO)CO  TO  666 

IF(I.EQ.J)EDERV=PI 

IF ( I .NE . J ) EDERV= DIMAG ( CUI*CTERM ) 

Z( I , J ) =BETA*  EDERV 
IF(ALPHA.EQ.O.DO)CO  TO  888 
CCCCCCCCCCCCC 

666  CT1= ( CZ ICAR-CZ J ) *CTERM 

CT2  =  EDELJ*  !  1 .  DO+CDLOG  (  CK/  (  CZ  ICAR-CZJ  r'l )  )  ) 

Z(I,J)=Z(I,J) +ALPHA*DREAL( CT1+CT2 ) 

C  IF(LAPOIS.EQ.l.AND.BETA.EQ.O.DO) Z(I,J)=10.0*Z( I,J) 

888  CONTINUE 

999  CONTINUE 

1000  CONTINUE 
RETURN 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  FIELD( RI ,N,XV1 ,XV2 , YV1 ,  YV2 ,  NM) 

IMPLICIT  COMPLEX* 16  (C) 

IMPLICIT  REAL*8(A-B,E-H,P-Z) 

DIMENSION  RI ( N ) ,XVl(NM) ,XV2(NM) ,YV1(NM) ,YV2(NM) 

COMMON/TYPE/LAPOIS , EPSR 

TPI=8.D0*DATAN(1.D0) 

IF( LAPOIS .EO. 1 )TPI=TPI*EPSR*8 .854D-12 
Cl*(  1 . DO  »0 . DO ) 

CK*(100.DO,O.DO) 

C 

C  READ  K=NO.  OF  POINTS  AT  WHICH  FIELD  AND  POTENTIAL  IS  TO  BE  COMPUTED 
C 

READ(92,*)XIN,YIN,XFIN,YFINf  K 

IF(K.E0.1)IJKL=1 

IF( K . EO • 1 ) GO  TO  78 

XDEL* ( XFIN-XIH ) /FLOAT ( K— 1 ) 

YDEL=( YFIN-YIN ) /FLOAT( K- 1 ) 

78  CONTINUE 

DO  50  J=1,K 

X= ( XIN+FLOAT (J-l)*XDEL)*l.D-06 
Y= ( YIN+FLOAT ( J  —  1 ) * YDEL ) *  1 . D-06 
IF  (  K  .  EO  •  1 .AND.J.NE. 1 ) CO  TO  50 
CZK=DCMPLX(X,Y) 

SUMP=0 . DO 
SUMX=0 . DO 
SUM Y=0 . DO 
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do  45  1=1, r: 

CZI=DCMPLX( XVI ( I )  , YVl ( I ) ) 

CZIP1=DCMPLX(XV2( I) ,YV2( I) ) 

EDELI=CDARS (CZIPl-CZI ) 

CUI  =  < CZIP1-CZX  J/EDELI 
CAP.C=  (CZK-CZIP1)/(CZK-CZ1) 

CTERMO=CDLOG(CARG)/CUI 
CTERM=  (CZK-CZI)  'CTERMO 

CTERM2=EDELI * ( Cl+CDLOC( CK/ ( CZK-C2 IF 1  , ) ) 

CWI=CTERM+CTERM2 

SUMP  =  SUMP+RI(  I)*DREAL(CVII) 

SUMX=SUMX+RI ( I ) *DR£AL( CTERMO ) 

SUMY=SUMY-RI ( I ) *  DIMAG ( CTERMO ) 

SUMF=SQRT<SUMX*SUMX+SUMY*SUMY) 

45  CONTINUE 

SUM P= SUMP/TP  I 
SUMX=SUMX/TPI 
SUMY=SUMY/TPI 
SUMF=SUMF/TPI 

C  WRITE (  18 ,55)CZK, SUMP, SUMX, SUMY, SUMF 

IF( LAPOIS . EO • 0 ) CO  TO  50 
CALL  POT EM ( X , Y , SP ) 

SUMP=SUMP+SP/TPI 
WRITE ( 18 , 49 ) X  ,  Y , SUMP 

49  FORMAT ( 3X , ' X= ' ,E11.5,2X,'Y='  , Ell . 5 , 3X , 1  TOTAL  POT= ’ , Ell .4 ,/) 

50  CONTINUE 

55  FORMAT ( IX, *  Z  =  * , 2E11 . 5 , 3X , ’ POT= ' , E10 . 4 , 3X , ’ EX= ’ ,E10.4,2X, 

6  'EY=‘ ,E10.4,3X, 'ETOT=' ,£10.4) 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  ELSV (A,B,C,N,DE,EP) 

IMPLICIT  REAL*8(A-H,P-Z) 

DIMENSION  A(N,N) ,B(N) ,C(N) 

DO  11  1=1, N 
B( I ) =0 . DO 
C( I)=0.D0 
DO  12  J=1,N 

12  C( I)=C( I)+A< I,J) 

11  A ( I , I ) =A ( I , I ) - 1 . DO 

DO  13  K=1 ,N 
DO  14  J=1 ,N 
B(J)=A<K,J) 

14  A(K,J)=0.D0 
A(  K  ,K )  *1 .  DO 
W=B( K ) +1 . DO 

IF{ ABS(W) .LT.EPJGO  TO  17 
DO  13  1=1, N 
Y=A( I ,K)/W 
DO  13  J=1 , N 

13  A(I,J)=A(I,J)-B(J)*Y 
DE=0 . DO 

DO  -15  J=1  ,M 
B ( J ) =0 . DO 
DO  16  1=1, M 

16  B(J)=B(J)+A(I,J) 

15  DE=CE+C(J)*B( J) 

RETURN 

17  DE=- 1 . DO 
RETURN 
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subroutine  pcten(x, y,pct) 
IMPLICIT  COMPLEX* 16  (C) 
IMPLICIT  PEAL*8  ( A-S , E-H , P- Z  ) 


C  THIS  PROGRAM  GIVES  THE  POTENTIAL  AND  THE 
C  GRADIENT  OF  THE  POTENTIAL  PRODUCED  BE  A 
C  CERTAIN  TWO  DIMENSIONAL  CHARGE  DISTRIBUTION 
C  AT  A  POINT  (X, Y) . 

C  NOTE  THAT  AS  THE  CHARGE  DISTRIBUTION  IS  CHANCED 
C  THIS  PP.OCRAM  SHOULD  BE  CHANCED  ACCORDINGLY. 

C 

C  THIS  PROGRAM  WILL  BE  CALLED  BY  THE  PROGRAM  NAMED 
C  ' DNEWPCISSCN ’  ONLY  IF  THE  PARAMETER  LAPOIS  IS  1  IN 
C  THAT  PROGRAM. 

C 

RHC=-3200 . DO 
AK=100 . DO 
Xl=-0. S0D-06 
Yl=-0 .04D-06 
X2- 0 . 5D-C6 
Y2  =  0 .04D-06 
CCCCCCCCCCCCCCCCCCC 

POT=(Y2-Y1)*(X2-X1)*DLOG(AK) 

CCCCCCCCCCCCCCCCCCC 

EDELl =DABS ( X2-X1 ) 

U1=(X2-X1 l/EDELl 
CCCCCCCCCCCCCCCCCCC 
FI1  =  ( Y2-Y1 ) 

POT=POT+£DELl *  FI I 
CCCCCCCCCCCCCCCCCCC 

CZ3=DCMPLX(X2, YI ) 

CZ4=DCMPLX(X2,Y2) 

CZ=DCMPLX( X , Y ) 

EDEL2=CDABS ( CZ4-CZ3 ) 

CU2=(CZ4-CZ3)/EDEL2 

CT2=( CZ-CZ3 ) *CDLQG{ ( CZ-CZ4 )/( CZ-CZ3 ) )/CU2 
CT2=CT2  +  EDEL2* ( 1 ,DO-CDLOG(  CZ-CZ4 ) ) 
CCCCCCCCCCCCCCCCCCCCC 
FI2=DREAL( CT2 ) 

FI31=(X1-X)*FI2 
POT=POT+EDELl*FI 2+FI 3 1/U 1 
CCCCCCCCCCCCCCCCCCCC 
XS  =  X2 

CALL  INTG (XS ,Y1 ,Y2,X,Y,RES32) 
CCCCCCCCCCCCCCCCCCCCCC 
FI32=RES32 
POT=POT  +  FI3  2/U1 
CCCCCCCCCCCCCCCCCCCC 

CZ5=GCMPLX(X1 , Yl ) 

CZ6=DCMPLX(X1 ,Y2) 

EDEL3=CDABS ( CZ6-CZ5 ) 

CU3  =  (CZ6-CZ5)/EDEL3 

CT4  1  =  (CZ-CZ5)*  CDLOC ( ( CZ-CZ6 )/ ( CZ-CZ5 ) )/CU3 
CT4l=CT41+EDEL3*( 1 . DO-CDLOC ( CZ -CZ 6 ) ) 
CCCCCCCCCCCCCCCCCCCCCCCC 

FI4l=(Xl-X)*DREAL( CT4 1 ) 

POT=POT-FI 4 1/U1 
CCCCCCCCCCCCCCCCCCCCCCC 


XS  =  X1 

CALL  INTC ( XS , Yl , Y2 , X , Y , RES4  2 ) 
CCCCCCCCCCCCCCCCCCCCCCCCCCCC 
FI42  =  P.ES42 
POT = POT- FI  4  2/Ul 
POT=RHO*POT 
RETURN 
END 

SUBROUTINE  INTCfXS , Y1 , Y2 ,  X  ,  Y  , RES ) 

IMPLICIT  REAL* 8 ( A-3 , E-H , P-Z ) 

PI=4 . D0*DATAN ( 1 . DO ) 

TP  =  2 . DO  *PI 
PO  2  =  PI/2 . DO 
P04=PI/4 . DO 
IF ( X-XS ) 20 , 10 , 15 

10  IF( Y.GE . Y2)RES=P04*( ( Y2- Y ) * ( Y2- Y ) - ( Yl- Y ) * ( Yl-Y) ) 

I F ( Y.LE. Y1 )RES=-P04* ( ( Y2-Y ) * { Y2-Y ) - ( Yl-Y) * ( Yl-Y) ) 

IF( Y.CT. Yl .AND. Y.LT. Y2 ) IM=1 

I F(  IM.EC.l  )P-ES=-P04*(  (  Y2-Y)  *  (  Y2-Y )  +  (  Yl-Y )  *  (  Yl-Y)  ) 

RETURN 

15  T1=(X-XS)*( Y2-Y15 

T2=( Y-Y2)*( Y-Y2)+{ X-XS)* (X-XS) 

T3= ( Y- Yl ) * ( Y- Yl )+ ( X-XS ) * ( X-XS ) 

T4= ( Y2- Y ) / ( X-XS ) 

T5=(Y1-Y)/(X-XS) 

IF(Y.GE.Y2)RES=(T 1-T2*DATAN ( T4 ) +T3*  DATAM ( T5 ) )/2.D0 
IF(Y.LE.Yl) RES  =  (T1-T2*  DATAN ( T4 ) +T3  *DATAN ( T5 )  )/2.D0 
IF( Y.CT. Yl. AND. Y.LT. Y2 ) IM=1 

I F ( IM.EO. 1 ) RES= ( T1 +T3* DATAN ( T5 ) -T2* DATAN ( T4 ) )/2 .DO 
RETURN 

20  T1=(XS-X)*( Y2-Y1) 

T2= ( Y2-Y ) * ( Y2-Y)  +  ( XS-X ) * ( XS-X) 

T3= ( Y-Yl ) * ( Y- Yl )+( XS-X ) * ( XS-X ) 

T4  =  ( Y2- Y)/ ( XS-X ) 

T5= ( Yl-Y)/ (XS-X) 

TERM=P02* ( (Y2-Y)*(Y2-Y)-( Yl-Y)* (Yl-Y) ) 

I F ( Y . CE . Y2 ) PES=TERM- ( T1-T2  *  DATAN ( T4 ) +T3  *DATAN ( T5 ) ) /2 . DO 
IF  (  Y.LE.  Yl  )P.ES  =  -TERM+(T2*  DATAN  (T4  ) -T1-T3*  DATAN  (  T5  )  )/2  .DO 
IF (Y.LT. Y2. AND. Y.CT. Yl ) IM=1 

I F  (  IM.EC.  1  )  TM  =  -PC2  ■’  (  (  Yl-Y)*  (Yl-Y)  +  (  Y2-Y)  *(  Y2-Y)  ) 

IF ( IM.EC.l ) RES=TM- ( T 1 +T3  *  DATAN ( T5 ) -T2* DATAN ( T4 ) J/2.D0 

PETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  GRAD ( X , Y , CX ,GY ) 
implicit  ccmplex*16  (c) 

IMPLICIT  REAL*8(A-B,E-K,P-Z) 

FHO=-3200 . DO 
C  2  =  DC  M  P  LX  (  X  ,  Y  ) 

XI =-0 . 50D-06 

Yl =-0 . 04D-06 

X2=0 .5D-06 

Y2  =  0 . 0  4  D-06 

CZ 1 =DCMPLX ( XI ,  Yl ) 

C2 2  =  DC."PLX  (X2  ,  Yl  ) 

CZ3=DCMPLX( XI , Y2) 

C  2  4  =  DCM P LX  (  X 2  ,  Y 2  ) 

EDEL1 =CrABS ( C2  2-C2 1 ) 


SAMPLE  INPUT/OUTPUT  FILE : 

The  following  is  the  input/output  file  for  the  long 
channel  MOSFET  problem  considered  in  Fig.  4.  The  results 
presented  in  the  following  are  plotted  in  Fig.  6  and  7. 
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'  y  uj  i  m  1  m  u  ■  mw  i  ■  |w^< 


j'V'  -. 

r-V-' 

.  V  .  ^ 

~  .  .  !  b  x 

*  ,*-  . * 

9.3  - C . 106  0.5 

0.106  12 

V\ 

1.0  0.0  5.3 

»  ,  '•  . ' 

-0.5  0.106  -C.5 

-0.105  12 

i 

l.G  0.0  0.3 
-0.5  -0.106  0.5 
0.3  1.0  0.0 

-0.106  25 

0.5  0.106  -0.5 

0.136  25 

1.0  3.22E-C6  2. 

0 

■_  ’ 

• 

■  ■  • 

-0.50  0.1059999 

0.50  0.105999  21 

|]3MB 

X= - . 50000E-06 

Y  =  0 . 10600E-06 

TOTAL 

PCT=  0 . 9552E+-00 

y[\ 

X=- . 45Q00E-06 

Y=C. 10600E-0& 

rrr'rr  \  r 

?OT  =  0.1016E-01 

'--"V 

i  .  V  . 

X=- . 4000CE-06 

Y=0. 10600E-06 

TOTAL 

?CT=  0 . 1005E+01 

p\  - 

& 

X=-. 35000E-06 

Y=0 . 10600E-06 

TOTAL 

POT  =  0.9841E-HQ0 

.. 

V  ''.  -  ‘ 

X=- . 30000E-06 

Y=0 . 1060CE-06 

TOTAL 

POT=  0 . 9704EtQ0 

V-V-' 

\  *  «  ' 

X=- . 25000E-06 

Y=0. 10600E-06 

TOTAL 

POT=  0 . 9691E+00 

$ 

X=-. 20000E-06 

Y=0. 10600E-06 

TOTAL 

POT-  0 . 9327E+00 

y  ■: 

X  =  -. I5000E-06 

Y-0 . 10600E-06 

TOTAL 

P9T=  0 . 1012E-I-01 

*  ' 

■/  .  ‘ 

X=- . 1Q000E-06 

Y=0 . 10600E-06 

TOTAL 

POT=  0 . 1059E-r01 

*  '* 

« 

X=- . 50000E-07 

Y=0. 10600E-06 

TOTAL 

POT=  0 . 1126E+01 

X=0 . OOOOQE+OO 

Y=0. 10600E-06 

TOTAL 

PCT=  0. 1216E+-01 

X=0 . 50000E-07 

Y=0 . 13600E-06 

TOTAL 

P0T=  0 . 1332E+01 

X=0 . 10000E-06 

Y=0. 10600E-06 

TOTAL 

PCT=  0 . 1478E+01 

a 

V  . 

X=0 . 15000E-06 

Y=0. 10600E-06 

TOTAL 

POT=  0.1661E-r01 

rH  .- 

r.  \ 

X=0. 20000E-06 

Y  =  0 . 10600E-06 

TOTAL 

POT=  0.1891E-r01 

<\\\ 

f.v. 

X-0. 25000E-06 

Y=0 . 10600E-06 

TOTAL 

POT=  0 . 2178E+01 

f.f? 

X-0 . 30000E-06 

Y-0 . 10600E-06 

TOTAL 

POT=  0 . 2534E+-0 1 

X=0 . 35000E-0& 

Y  -0 . 10600E-06 

TOTAL 

POT=  0 . 2936E-0 1 

X=0 . 40000E-06 

Y  -  0 . 1 060 0E- 06 

TOTAL 

?0T=  0.3575E+01 

X  =  0 . 45  300E-06 

Y=C. 1060CE-06 

TOTAL 

P0T=  0 .  133~Ei-01 

* 

X=0. 500CCE-06 

Y=0 . ICbOCE-OS 

TOTAL 

POT=  0.5556E-91 

•3 


I  -  W.  m  A  «.  *  m  f 


y^.^\V‘j 


o  m  o 


.1  0.0  3.3 

2.3  0.106  - 2 . 5  -C.106  12 
.0  0.0  0.9 

2.5  -0.106  2.5  -0.106  50 

1.0  0.0 

0.106  -2.5  0.106  50 

0.22E-06  2.0 

2.50  0.1059999  2.50  0.105999  21 


A  —  — . 250 COE  —  05 

Y-0 . 10600E-06 

TOTAL 

POT -  0.9233ErQ0 

X- - . 2250CE-05 

Y  =  0 . 10600E-06 

TOTAL 

?CT=  0 . 7867E-0  0 

X^-. 20000E-05 

Y=0 . 106C0E-06 

TOTAL 

POT=  0 . 6475E-00 

X=- . 17500E-05 

Y=C . 10600E-06 

TOTAL 

POT=  0.5955E-r00 

X=  - . 15000E-05 

Y=0 . 1C600E-06 

TOTAL 

POT  =  0 . 5764St-00 

X=  - . 12500E-05 

Y-0 . 1060CE-06 

TOTAL 

POT=  0 . 5694E+-00 

X-  - . 10000E-05 

Y=0 . 10600E-06 

TOTAL 

POT=  0 . 5669E+00 

X-- . 75000E-06 

Y=0 . 10600E-06 

TOTAL 

POT=  0 . 5659E-00 

X=  - . 50CC  0E-O6 

Y=0. 10600E-06 

TOTAL 

POT=  0 . 5o56E-i-00 

X=  -  . 25000E- 06 

Y=0 . 10600E-06 

TOTAL 

POT=  0.5655Ei-00 

X=C .  OCOOOEt-OO 

Y=0 . 10600E-06 

TOTAL 

POT=  0.5656E+C0 

X=0 -2500CE-06 

Y=0. 10600E-06 

TOTAL 

POT=  0 . 5659Et00 

X=0. 50000E-06 

Y=0. 10600E-06 

TOTAL 

?0T=  0.5667E+00 

X=0 . 75000E-06 

Y-0 . 10600E-06 

TOTAL 

PCT=  0 . 5692E+00 

X=0. iCCOOE-05 

Y=0 . 10600E-06 

TOTAL 

POT=  0 . 5753E+-00 

X=0. 125C0E-05 

Y=0. 10600E-06 

TOTAL 

POT=  0 . 5940E->-00 

X=0. 15000E-05 

Y=0 . 10600E-06 

TOTAL 

POT=  0.6434E-f-00 

X--0  .  175C0E-C5 

Y-0 . 10600E-06 

TOTAL 

POT=  0 . 773SE+00 

X=0 . 200C0E-05 

Y=0 . 1060CE-06 

TOTAL 

POT=  0.1147Eb01 

X=0. 22500E-05 

Y -  0 . 1060CE-06 

TOTAL 

POT -  0 . 2158E-01 

X-  0 . 25  0OOE-O5 

Y--0 .10600E-06 

TOTAL 

POT  =  0 . 5439E-01 
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Fig.  2.  The  potential  in  R  is  produced  by  the  impressed 

volume  charges  p  and  the  equivalent  surface  charges  a. 

v  + 

The  surface  charges  are  on  C  ,  (just  outside  of  bounding 
surface  C ) . 
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Surface  Potential  Variation  Along  the  Channel  for 
1  n  m  Device 


Threshold  Voltage  Versus  Channel 


